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A  new  approach  to  numerical  simulation  of  liquid  water  distribution  in  channels  and  porous  media  includ¬ 
ing  gas  diffusion  layers  (GDLs),  catalyst  layers,  and  the  membrane  of  a  proton  exchange  membrane  fuel 
cell  (PEMFC)  was  introduced  in  this  study.  The  three-dimensional,  PEMFC  model  with  detailed  thermo¬ 
electrochemistry,  multi-species,  and  two-phase  interactions.  Explicit  gas-liquid  interface  tracking  was 
performed  by  using  Computational  Fluid  Dynamics  (CFD)  software  package  FLUENT®  v6.2,  with  its  User- 
Defined  Functions  (UDF)  combined  with  volume-of-fluid  (VOF)  algorithm.  The  liquid  water  transport  on 
a  PEMFC  with  interdigitated  design  was  investigated.  The  behavior  of  liquid  water  was  understood  by 
presenting  the  motion  of  liquid  water  droplet  in  the  channels  and  the  porous  media  at  different  time 
instants.  The  numerical  results  show  that  removal  of  liquid  water  strongly  depends  on  the  magnitude  of 
the  flow  field.  Due  to  the  blockage  of  liquid  water,  the  gas  flow  is  unevenly  distributed,  the  high  pressure 
regions  takes  place  at  the  locations  where  water  liquid  appears.  In  addition,  mass  transport  of  the  species 
and  the  current  density  distribution  is  significantly  degraded  by  the  presence  of  liquid  water. 

©  2009  Elsevier  B.V.  All  rights  reserved. 


1.  Introduction 

Actual  PEMFC  operation  is  a  complex  process  including  trans¬ 
port  of  mass,  momentum,  energy,  species  and  charges  that  take 
place  simultaneously.  Furthermore,  a  PEMFC  is  a  multi-part  device 
that  comprises  of  collectors,  flow  channels,  GDLs,  catalyst  layers 
and  the  PEM  membrane.  During  a  PEMFC  operation,  hydrogen  is 
consumed  in  the  anode;  oxygen  is  consumed  and  water  is  pro¬ 
duced  in  the  cathode;  and  a  current  is  generated  from  anode  side 
to  cathode  side  via  electric  load.  Water  product  exists  in  the  forms 
of  water  vapor  and/or  liquid  water  depending  on  physical  condi¬ 
tions  such  as  pressure,  temperature,  etc.,  leading  to  a  gas-liquid 
flow  in  the  porous  media  and  the  channels.  Such  water  signifi¬ 
cantly  influences  the  fuel  cell  operation  and  performance  as  well: 
water  content  increases  ionic  conductivity  of  the  membrane,  lead¬ 
ing  to  a  high  current  density  in  a  specific  volume  of  the  cell;  liquid 
water  may  block  the  gas  transport  inside  the  channel  and  porous 
media  and  cause  flooding  phenomenon  and  therefore,  decreasing 
the  cell  performance.  Therefore,  water  management,  especially  liq¬ 
uid  water  management,  to  which  many  engineers  and  scientists 
have  recently  paid  particular  attention,  has  been  a  critical  challenge 
for  high-performance  fuel  cell  design  and  optimization. 

In  order  to  fully  understand  physical  and  electrochemical  pro¬ 
cesses  occurring  in  a  PEMFC,  and  especially  water  management, 


*  Corresponding  author.  Tel.:  +1  519  253  3000x2630;  fax:  +1  519  973  7007. 
E-mail  address:  bzhou@uwindsor.ca  (B.  Zhou). 


numerical  modeling  is  a  strong  method  that  has  been  much  con¬ 
sidered  in  recent  years.  As  mentioned  above,  PEMFC  operation 
involves  simultaneous,  complex  processes  and  hence,  the  devel¬ 
opment  of  3-D  models  for  a  complete  PEMFC  that  considers  all 
parts,  all  components  and  all  phenomena  is  a  big  challenge.  For¬ 
tunately,  with  the  developments  of  high  performance  computing 
and  advanced  numerical  algorithms,  it  has  allowed  researchers 
to  model  PEMFC  systems  as  well  as  individual  components  with 
greater  fidelity  than  ever  before.  Recently,  three-dimensional  Com¬ 
putational  Fluid  Dynamics  (CFD)  models  have  been  developed  by 
taking  full  advantage  of  different  commercial  CFD  software  pack¬ 
ages  such  as  Fluent®,  CFX®,  Star-CD®,  and  CFDRC®,  etc.  Typical 
models  that  simultaneously  considered  the  electrochemical  kinet¬ 
ics,  current  distributions,  fluid  dynamics,  and  multi-component 
transport  accounted  for  single-phase  and  two-phase  were  con¬ 
ducted  by  Wang  and  co-workers  [1-4].  In  these  studies,  liquid  water 
saturation  was  addressed  to  characterize  the  two-phase  flow  in 
the  cathode  channel  and  cathode  GDL.  The  CFD  software  Fluent 
was  also  used  to  simulate  fluid  flow  in  PEMFC  channel  by  Dutta 
et  al.  [5]  or  predict  liquid  water  saturation  with  experimentally 
measured  capillary  functions  by  Ye  and  Nguyen  [6].  Djilali  et  al.  pre¬ 
sented  a  multi-phase,  multi-component  3-D  model  with  heat  and 
mass  transfer,  where  liquid  water  transport  in  GDLs  was  numeri¬ 
cally  presented  [7,8].  Another  model  to  treat  liquid  water  formation 
and  transport  in  a  PEMFC  was  introduced  by  Mazumder  and  Cole 

[9] .  Other  approach  to  deal  with  two-phase  flow  in  GDLs  based 
on  thermodynamic  equilibrium  conditions  was  given  by  Vynnycky 

[10] .  By  using  this  approach,  the  location  of  the  interface  between 
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Nomenclature 

a 

water  activity,  dimensionless 

^surf 

reactive  surface  area,  m2 

heat  transfer  surface  area,  m2 

cp 

specific  heat  capacity,  J  kg-1 K-1 

Q 

species  concentration  z,  kmol  m-3 

Di 

diffusion  coefficient  of  species  i  in  gas  mixture, 

m2  s-1 

F 

Faraday  constant,  9.6487  x  107  Ckmol-1 

h 

convective  heat  transfer,  W  m-2 I<-1 

^h2o 

the  enthalpy  of  formation  of  water  vapor,  N  m  kg-1 

/ 

current  density,  Am-2 

Gve 

average  current  density,  A  m-2 

J 

mass  flux,  kgm-2  s-1 

/<eff 

effective  thermal  conductivity,  Wm-1 I<-1 

Mi 

molecular  weight  of  species  i  in  gas  mixture, 

kg  kmol-1 

nd 

electro-osmotic  drag  coefficient,  dimensionless 

Hf 

charge  number  of  the  sulfonic  acid  ion 

P 

pressure,  Pa 

Pi 

partial  pressure  of  species  i,  Pa 

R 

universal  gas  constant,  8314J  kmol-1  K-1 

Q 

heat  rate,  W 

ftcat,  Ran 

volumetric  current  density,  Am-3 

S 

source  term 

t 

time,  s 

T 

temperature,  I< 

U ,  V ,  W 

velocities  in  x,  y ,  and  z  directions,  respectively,  m  s-1 

Voc 

open-circuit  potential,  V 

Vcell 

cell  potential,  V 

Vref 

reference  potential,  V 

V 

volume,  m3 

X, 

mole  fraction  of  species  i,  dimensionless 

Vi 

mass  fraction  of  species  i,  dimensionless 

Greek  symbols 

a 

transfer  coefficient, 

P 

the  factor  accounts  for  energy  release 

£ 

porosity 

0 

phase  potential,  V 

y 

concentration  dependence 

Kp»  Yt 

exponent  factors 

T) 

overpotential,  V 

(p 

relative  water  content 

K 

surface  curvature 

H 

gaseous  permeability,  m2 

liquid  permeability,  m2 

Tp 

hydraulic  permeability,  m2 

~t(p 

electrokinetic  permeability,  m2 

X 

water  content 

jl 

dynamic  viscosity,  kgm-2  s-1 

P 

density  of  gas  mixture,  kg  m-3 

Pi 

density  of  species  i,  kg  m-3 

o 

phase  conductivity,  £2-1  m-1 

V 

reaction  rate,  kmol  m2  s-1 

00 

excess  coefficient 

X 

surface  tension  coefficient,  N  m-1 

0w 

contact  angle, 0 

Subscripts  and  superscripts 

an 

anode 

cat 

cathode 

e 

electrochemical  reaction 

eff 

effective 

g 

gas  phase 

h2 

hydrogen 

h2o 

water 

i 

species  i 

in 

inlet 

1 

liquid  phase 

m 

membrane  phase 

out 

outlet 

n2 

nitrogen 

02 

oxygen 

ref 

reference 

s 

solid  phase 

sat 

saturated 

surf 

surface 

w 

water  vapor 

one  phase  and  the  other  phase  was  discovered  from  its  numerical 
results. 

Two-phase  flow  models  that  account  for  liquid  water  forma¬ 
tion  in  the  porous  media  including  GDLs  and  catalyst  layers  were 
developed  in  several  ways.  A  typical  approach  using  the  network 
model  to  investigate  the  effective  diffusivity  and  water  saturation 
was  developed  early  by  Nam  and  Kaviany  [11].  In  this  study,  the 
effects  of  the  fiber  diameter,  porosity,  and  capillary  pressure  on  the 
water  saturation  and  the  cell  performance  are  also  presented.  The 
effective  permeability  and  capillary  pressure  in  porous  GDLs  were 
also  addressed  by  using  capillary  pore  network  simulations  given 
by  Djilali  et  al.  [12].  Consequently,  the  pore  network  method  has 
been  used  to  address  the  liquid  water  behavior  in  a  PEMFC  by  Fowler 
and  co-workers  [13],  Prat  and  co-workers  [14],  and  Pasaogullari  and 
Wang  [15]. 

A  new,  comprehensive  approach  used  to  address  the  behavior 
of  liquid  water  transport  inside  PEMFCs  was  early  developed  by 
Zhou  and  co-workers  [16-18]  that  considers  volume-of-fluid  (VOF) 
interface  tracking  between  liquid  water  and  gas.  This  study  also  con¬ 
ducted  two  more  studies  that  dealt  with  liquid  water  in  serpentine 
and  straight  parallel  fuel  cell  stacks.  The  results  showed  that  dif¬ 
ferent  designs  of  gas  diffusion  layers  (GDLs)  and  flow  channels  will 
affect  the  liquid  water  flow  patterns  significantly,  thus  influencing 
the  performance  of  PEMFCs.  In  recent  years,  Djilali  and  co-workers 
[19]  presented  a  2-D,  numerical  investigation  of  the  dynamic  behav¬ 
ior  of  liquid  water  entering  a  PEMFC  channel  through  a  small  hole 
that  was  assumed  as  a  GDL  pore.  The  3-D  simulation  models  focus¬ 
ing  on  characteristics  of  water  droplets  and  water  film  motion  in 
the  flow  channels  and/or  GDL  were  also  developed  by  using  VOF 
method  by  Djilali  and  co-workers  [20,21  ]  Zhan  et  al.  [22],  Quan  and 
Lai  [23],  Golpaygan  and  Ashgriz  [24]. 

In  these  studies,  the  water  behavior  can  be  investigated,  giving 
very  useful  insights  about  water  management  in  the  channels,  how¬ 
ever  the  electrochemical  reactions  and  the  effects  of  liquid  water 
on  mass  and  heat  transport  that  take  place  in  an  actual  PEMFC  have 
not  been  considered. 

Therefore,  a  general  3-D  PEMFC  model  including  the  entire 
serpentine  flow  channels,  GDLs,  catalyst  layers,  membrane  and  col¬ 
lectors  all  with  detailed  physics  included  (e.g.,  multi-phase  with 
VOF  interface-tracking  between  gas-phase  and  liquid-water-phase, 
multi-components,  heat  and  mass  transfer,  electrochemical  reac¬ 
tions,  and  water-phase-change  effects)  has  been  developed  by  Le 
and  Zhou  [25].  In  this  study,  as  a  pioneering  work  in  this  research 
area,  the  mathematical  model  is  applied  to  a  PEMFC  with  interdigi- 
tated  flow  channels.  A  unified  method  is  used  to  treat  the  transport 
of  water  in  different  zones:  the  species  transport  equation  of  water 
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Fig.  1.  (a)  Schematic  diagram  of  the  PEMFC  and  (b)  geometry  and  computational 
mesh  of  the  interdigitated  channel. 


vapor  coupled  with  the  water  transport  through  the  membrane 
(including  the  electro-osmotic  drag  and  the  back  diffusion  effects) 
is  used  for  different  zones  with  different  source  terms.  Further¬ 
more,  we  also  use  VOF  method  to  track  liquid  water  and  treat  water 
vapor  as  one  of  the  gaseous  species.  By  this  way,  liquid  water  motion 
and  behavior  inside  the  porous  media  are  implicitly  investigated  in 
the  flow  interdigitated  channel  with  baffles.  The  simulation  results 
would  provide  an  in-depth  understanding  of  liquid  water  trans¬ 
port  and  its  effects  on  fuel  cell  performance  in  an  interdigitated 
flow-field  PEMFC. 


2.  Mathematical  model 

Fig.  1  shows  schematic  structure  of  the  3-D  PEMFC  used  in  this 
study.  The  mathematical  model  considers  all  the  parts  shown  in 
Fig.  la,  the  computation  grid  of  interdigitated  channel  is  shown  in 
part  b,  with  geometrical  parameters  listed  in  Table  1.  The  compu¬ 
tation  domain  consists  of  current  collectors,  flow  channels,  GDLs, 
catalyst  layers  and  the  membrane.  The  channels  are  interdigitated. 

2.1  Model  assumptions 

The  assumptions  used  in  developing  the  model  are: 

1.  Ideal  gas  law  was  employed  for  gaseous  species. 

2.  The  fluid  flow  in  the  fuel  cell  was  laminar  due  to  the  low  flow 
velocities  and  the  small  size  of  gas  flow  channels. 

3.  The  porous  media  including  membrane,  catalyst  layers  and  GDLs 
were  considered  to  be  isotropic. 

4.  The  fuel  cell  cooling  was  controlled  by  forced  convection  heat 
transfer. 


Table  1 

Geometrical  properties  and  operation  conditions. 


Physical  properties  or  parameters 

Value 

Channel  width 

0.001  m 

Channel  height 

0.001  m 

Membrane  thickness 

50  x 10-6  m 

GDL  thickness 

300  x  10“6  m 

Catalyst  layer  thickness 

10  x 10-6  m 

Anode  and  Cathode  Inlet  Temperature 

300  K 

Anode  inlet  excess  coefficient,  coa 

3 

Cathode  inlet  excess  coefficient,  coc 

3 

Open-circuit  voltage,  V0c 

1.15  V 

Gas  Constant,  R 

8.314  Jkmol-1  K_1 

Anode  exchange  current  density/reference 

7  x  1010  Akmol-1 

hydrogen  concentration,  ^7(6^  ^ 

Cathode  exchange  current  density/reference 

7  x  105  Akmol-1 

oxygen  concentration,  ^|[/(CQe2f)ycat 

Anode  transfer  coefficient,  aan 

0.5 

Cathode  transfer  coefficient,  a^at 

0.5 

Anode  concentration  dependence,  yan 

0.5 

Cathode  concentration  dependence,  yCat 

1.0 

Factor  accounts  for  energy  release,  ft 

0.5 

Membrane  porosity,  emem 

0.5 

Diffusion  layer  porosity,  egd\ 

0.5 

Catalyst  layer  porosity,  eCataiyst 

0.5 

Permeability  of  porous  media,  r 

1.76  x  10-11  m2 

Contact  angle,  0W 

90° 

Surface  tension,  x 

0.065  Nm-1 

2.2.  Governing  equations 

In  this  study,  the  PEMFC  model  includes  the  principles  of  con¬ 
servation:  mass  and  momentum  for  solving  fluid  flow,  energy  for 
solving  heat  generation  and  transfer,  species  for  solving  species 
transport,  and  charge  for  solving  current  density  distribution.  In 
addition,  a  volume  fraction  equation  was  used  to  track  gas-liquid 
interfaces.  These  governing  equations  are  presented  in  the  follow¬ 
ing  forms: 


Continuity  and  momentum  equations 

^  +  V.(ept>)  =  Sm 

(i) 

9 

—{spv)  +  W{epvv)  =  -sVp  +  V[£/xVp]  +  Sv 
ot 

(2) 

The  source  terms  for  the  continuity  and  momentum  equation 
used  in  the  model  for  different  regions  of  the  fuel  cell  are  given 
in  Table  2.  For  the  different  fluid  zones,  different  source  terms  are 
used  to  describe  the  flow  of  the  fluid  through  a  porous  media  by 
using  viscous  loss— Darcy’s  drag  force  [26].  In  addition,  the  gravity 
and  surface  tension  forces  were  also  considered  in  the  momentum 
equation  source  term. 

Mixture  properties:  The  gas  and  liquid  phases  are  considered  as 
a  two-phase  mixture  flow.  For  the  mixture  flow  including  gas  and 
liquid,  the  volume  fractions  of  the  gas  and  liquid  phases  in  the 
computational  cell  are  introduced.  In  each  control  volume,  the  sum 
of  volume  fractions  of  two  phases  is  unity.  Then,  mixture  properties 
are  defined  as  in  Table  3. 

Volume  fraction  equation  of  liquid  water 

9 

gp(esipi)  +  V(esipiPi)  =  Ss  (3) 

The  tracking  of  the  interface  between  the  phases  was  accom¬ 
plished  by  the  solution  of  a  continuity  equation  for  the  volume 
fraction  of  one  (or  more)  of  the  phases.  The  motion  of  the  interface 
between  two  immiscible  liquids  (namely,  the  gas  and  liquid  water) 
of  different  density  and  viscosity  in  the  VOF  method  was  defined 
by  volume  fraction  of  liquid  water  S\  and  volume  fraction  of  the 
gaseous  phase  sg  [26].  In  this  model,  the  volume  fraction  equation 
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Table  2 

The  source  terms  of  governing  equations. 


Governing  equation 


Volumetric  source  terms  and  location  of 
application 


Conservation  of  mass 


Volume  fraction 

Conservation  of  momentum 


Conservation  of  energy 


Hydrogen  transport 
Oxygen  transport 

Water  vapor  transport 

Conservation  of  Charge 


For  gas  channels,  GDL  and  the  membrane 
Sm  =  0 


For  anode  catalyst  layer: 

C  _  Mh2  D  (  ndMH20 

Jm  —  2F~  ^an  —  1 

For  cathode  catalyst  layer: 


^2°  A  D 

F  J 


C  _  U2  D 

Jm  —  4f^^cat 

For  all  parts: 

Ss  =  Gv 


mH20  p  (  ndMH- 

H - 2 F~^cat  "h  I  - ; 


For  gas  channels: 
5V  =  0 


For  GDLs  and  void  of  catalyst  layers: 
Sv  =  PS-fs2J  +  XK(^|) 

For  membrane: 

Sv  =  PS-  + 

For  current  collectors: 

St  =  (/2/<ts) 

For  gas  flow  channel: 

St  =  rwhL 

For  GDL: 

ST  =  -4r  +  rw^L 
°s 

For  membrane: 

ST  =  4*  +  rwhL 


(4ff  +  4f) 


For  catalyst  layer: 

St  =  ^Fan.cat  +  J2  | 

For  anode  catalyst  layer: 

C  Mh2  D 

■JH2  = - 2F~Kan 

For  cathode  catalyst  layer: 

Mn 

S02  =-TR« 

For  anode  catalyst  layer: 

/  ndMH  o 

Sh20  =  -  I  — F - 

For  cathode  catalyst  layer: 

c  _  Mh2° 

^h20  —  2  F 

For  anode  catalyst  layer: 

Ses  =  — Fan 
Sem  —  Fan 

For  cathode  catalyst  layer: 
Ses  =  Fcat 
Sem  =  — Feat 
For  other  parts: 

Ses  —  0 
Sem  =  0 


+  rwhL 


)«a, 

st  lav 
Rcat+(^)l 


is  solved  for  liquid  water  phase.  The  volume  fraction  of  the  gas 
phase  is  automatically  computed  based  on  the  relative  equation 
listed  in  Table  3 
Energy  conservation  equation 


(pepwf  +  (pcp)eff(vVT)  =  V  I  ke[[VT  -  Y^hjjj  +  (r  ■  V)  1  +ST 


(4) 


right-hand  side  of  Eq.  (4)  represent  energy  transfer  due  to  conduc¬ 
tion  taking  place  in  solid  and  fluid  zones,  species  diffusion,  and 
viscous  dissipation,  respectively  [26].  For  the  collectors  what  are 
pure  solid  materials,  the  species  diffusion,  viscous  dissipation  and 
convection  term  are  negligible.  For  the  channel  and  porous  media, 
those  three  terms  were  considered  since  the  fluid  flow  is  presented 
in  such  regions.  Note  that  the  porosity  specified  in  the  energy  equa¬ 
tion  determines  different  zones,  depending  on  its  values:  1  for  the 
pure  fluid  zone  (the  channels),  0  for  pure  solid  zone  (the  collectors 
and  ribs),  and  between  0  and  1  for  porous  zones.  The  source  term 
in  the  energy  conservation  includes  heat  from  chemical  reactions 
(only  in  the  catalyst  layers),  Ohmic  heating  due  to  Ohmic  resistance 
of  solid  zones,  and  phase  change  (condensation/evaporation  pro¬ 
cesses).  Hence,  the  sources  terms  of  energy  equation  for  different 
layers  are  different  that  listed  in  detail  in  Table  2. 

Species  transport  equations 

The  species  transport  equations  are  generally  in  the  following 
form: 

^(spYi)  +  V  ■  (epvYj)  =  Dj  mV2(pYj)  +  S,-  (5) 

The  model  predicts  the  local  mass  fraction  of  each  species,  Yif 
through  the  solution  of  a  convection-diffusion  equation  for  the 
ith  species,  where  Di  m  is  the  diffusion  coefficient  for  species  i  in 
the  mixture  shown  in  Table  3.  The  source  term  in  the  transport 
equation  is  shown  in  Table  2. 

During  the  operation,  the  H+  ions  move  from  the  anode  to  the 
cathode  and  also  pull  water  molecules  with  their  movement.  This 
is  known  as  the  electro-osmotic  drag  effect.  Physically,  the  water 
transport  rate  through  the  membrane  from  anode  to  cathode  by 
electro-osmotic  drag  is  computed  as  a  function  of  volumetric  cur¬ 
rent  density  and  drag  coefficient  proposed  by  Springer  et  al.  [28] 
for  a  Nafion  membrane.  In  Springer’s  equation,  the  drag  coefficient 
is  defined  to  be  a  function  of  (the  water  content  inside  the  polymer 
membrane).  These  parameters  are  listed  in  Table  3. 

Conservation  of  charge 

The  current  transport  of  electrons  through  the  solid  phase  and 
ions  through  the  membrane  phase  was  represented  by  the  follow¬ 
ing  equations  [26,27]: 

V  •  (crs V0S)  =  Ses  (6) 

V  •  (crmV0m)  =  Sem  (7) 

The  volumetric  source  terms  Ses  and  Sem  are  defined  as  the  volu¬ 
metric  transfer  currents.  The  source  term  in  the  equation  of  charge 
is  shown  in  Table  2. 


The  source  terms  representing  the  transfer  currents  were  calcu¬ 
lated  by  using  Butler-Volmer  equation  [23]: 


(9) 


The  energy  balance  in  terms  of  the  temperature  change  was  also 
considered.  In  the  multi-phase  model,  the  energy  equation  is  also 
shared  among  the  phases  [26,27].  To  specify  parameters  in  the  mix¬ 
ture  including  gas  and  liquid  phases  in  porous  media,  the  effective 
properties  were  determined  (Table  3).  The  first  three  terms  on  the 


The  volumetric  transfer  current  R  is  driven  by  the  activation 
overpotential  ij,  which  is  the  potential  difference  between  solid  and 
membrane  phases  [26,27].  Note  that  the  reference  potential  of  the 
electrode  Vref  is  zero  on  the  anode  side  and  is  equal  to  the  open 
circuit  voltage  on  the  cathode  side.  The  cell  potential  is  then  the  dif¬ 
ference  between  the  cathode  and  anode  solid  phase  at  the  terminal 
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Table  3 

Physical,  species,  heat,  mass  and  charge  transport  properties. 
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Value/formulation 

Physical  parameter 

Mixture  density 

p  =  sipi  +  (l-si)pg 

Volume  fractions  relation 

Si+Sg  =  l 

Mixture  viscosity 

p  =  SifM  +  (\ -Si)pg 

Local  mass  average  velocity  of  mixture 

v  _  S\P\V\  +sgPgyg 
siPi  +sgPg 

Effective  volumetric  heat  capacity  of  porous  media 

(pCp)eff  =  epfCPif  +  (l  —  e)psCp,s 

Effective  thermal  conductivity  of  porous  media 

/<eff  =  £/<s  +  (l  -e)kf 

Heat  transport  parameter 

Heat  transfer  rate  due  to  forced  convection 

Qconvection  =  hAs(Ts  —  T00)(W) 

Heat  generation  rate  for  the  fuel  cell 

Species  transport  parameter 

Qgenerated  =  K  1.25  —  Vcell  )(  W) 

Water  transport  through  the  membrane  by  electro-osmotic  drag 

ndMH20  D 

1TIh20  —  F  Kcat 

Drag  coefficient 

nd  =  (2.5/22)A, 

Membrane  ionic  conductivity 

<xm  =  £(0.514A.  -  0.326)e1268((1/303)-(1/7')) 

Water  content 

Water  activity 

k  =  0.043  +  17.18a  -  39.85a2  +  36a3  (0 
k  =  14  +  1.4(a-  1)  (1 

k  =  16.8  (a 

a=  +sY 

1  sat  1  sat 

<  a  <  1) 

<  a  <  3) 

>3) 

Saturated  pressure 

Psat  =  101325 

^|q(— 2. 1794+0. 02953*Cr— 273.17)— 9. 1837*10-5*(T— 273. 17)2 +1. 4454*10-7*(T— 273. 17)3) 

Species  concentrations 

Q-Ytp/M, 

Mass  transport  parameter 

Phase  change  rate  (Condensation/Evaporation) 

rw  =  cr  max  [(1  -  s, ) M,h0,  -s,pi] 

Charge  transport  parameter 

Activation  overpotential 

r]  =  <Ps  ~  0m  —  Vref 

Cell  potential 

i^cell  =  0s, cat  _  0s, an 

Average  current  density 

/ave=W  Pand=l/  Rcatd 

J  ^an  J  Vcat 

collectors  (two  ends  of  the  cell  collectors  which  are  connected  to 
the  external  electric  circuit  from  both  electrodes). 

In  order  to  satisfy  the  conservation  of  charge,  the  total  current 
of  either  electrons  or  protons  coming  out  from  the  anode  catalyst 
layer  must  be  equal  to  the  total  current  coming  into  the  cathode 
catalyst  layer  and  must  be  equal  to  the  total  current  caused  by  the 
proton  movement  through  the  membrane  [28]: 

[  Rand  =  [  Rcatd  (10) 

3Van  Jv cat 

2.3.  Water  transport  and  its  effect  on  the  properties  of  porous 
media 

Water  is  formed  in  the  cathode  catalyst  layer  by  electrochemi¬ 
cal  reactions— an  amount  of  oxygen  is  consumed  and  an  amount  of 
water  is  produced.  Due  to  the  proton  movement  from  the  anode  to 
the  cathode  through  the  membrane,  water  molecules  are  pulled 
with  the  protons  by  a  force  called  electro-osmotic  drag  [26,29]. 
Additionally,  water  may  diffuse  through  the  membrane  due  to  the 
concentration  differences.  The  net  water  flux  through  the  mem¬ 
brane  results  in  a  water  balance  between  the  electro-osmotic  drag 
of  water  (from  anode  to  cathode)  and  back  diffusion  (from  cathode 
to  anode) 

Jw  =  JS?  +J^{  =  nd  -  Di  mV(fiYi)  ( 11 ) 


In  the  present  study,  the  membrane  water  diffusivity  and  mem¬ 
brane  ionic  conductivity  were  calculated  by  the  expressions  given 
by  Springer  et  al.  [28,29]  and  are  listed  in  Table  3.  It  is  almost  a 
standard  model  and  therefore  the  readers  can  check  the  original 
references  for  more  details. 

2.4.  Volume  of  fluid  (VOF)  model 

The  VOF  technique  was  implemented  in  the  channels  and  porous 
media  (including  GDLs  and  catalyst  layers)  [26].  Interface  between 
gas  and  liquid  (two-phase  flow)  is  tracked  by  the  volume  fraction  of 
liquid  water  in  the  computational  cell  volume.  In  the  VOF  approach, 
the  source  terms  of  continuity  and  momentum  equations  used  in 
the  porous  media  include  the  effects  of  surface  tension,  wall  adhe¬ 
sion  and  capillary  water  transport  phenomenon. 

2.4.1.  Geometric  reconstruction  scheme 

The  geometric  reconstruction  PLIC  scheme  (piecewise  linear 
interface  construction)  was  employed  because  of  its  accuracy  and 
applicability  for  general  unstructured  meshes,  compared  to  other 
methods  such  as  the  donor-acceptor,  Euler  explicit,  and  implicit 
schemes.  A  VOF  geometric  reconstruction  scheme  is  divided  into 
two  parts:  a  reconstruction  step  and  a  propagation  step.  Details 
can  be  found  in  [26,30]. 

2.4.2.  Implementation  of  surface  tension 

The  addition  of  surface  tension  to  the  VOF  method  is  modeled 
by  a  source  term  in  the  momentum  equation.  The  pressure  drop 
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across  the  surface  depends  upon  the  surface  tension  coefficient  x 
[26] 

Ap-xr,x(l  +  l)  (12) 

where  x  denotes  the  surface  tension  coefficient,  k  is  the  surface 
curvature;  R \  and  R2  are  the  two  radii,  in  orthogonal  directions, 
to  measure  the  surface  curvature.  The  surface  tension  force  can  be 
written  in  terms  of  the  pressure  jump  across  the  interface,  which  is 
expressed  as  a  volumetric  force  F  added  to  the  momentum  equation 


F  = 


pVsj 

xkU^)J2 


(13) 


where  the  term  (p/(pi  +  P2)/2)  is  used  to  improve  the  Continuum 
Surface  Force  (CSF)  method’s  capability  of  modeling  surface  ten¬ 
sion  in  the  case  for  fine  meshes  at  high  density  ratio  interfaces. 
This  method  is  called  “density-scaling  of  the  CFS”.  It  helps  to  main¬ 
tain  a  constant  interface  thickness  when  the  interface  between  a 
dense  phase  (namely,  liquid  water)  and  a  light  phase  (namely,  air) 
is  tracked  by  using  VOF  algorithm  [32].  The  source  terms  for  differ¬ 
ent  regions  of  the  fuel  cell  are  given  in  Table  2.  The  surface  curvature 
k  can  be  defined  in  terms  of  the  divergence  of  the  normal  unit  vector 
of  the  interface  h 


/c  —  V  h  —  V  •  (hssCOS 0$s  +  tssSin 0SS)  (14) 

where  h  is  the  unit  vector  normal  to  the  interface  between  two 
phases  near  the  solid  surfaces,  hss  is  the  unit  vector  normal  to  the 
solid  surfaces,  fss  is  the  unit  vector  tangential  to  the  solid  surfaces, 
and  0SS  is  the  static  contact  angle  at  the  solid  surfaces.  For  the  elec¬ 
trode  surfaces  with  different  wetabilities,  different  static  contact 
angles  could  be  assigned,  and  different  contact  angles  could  result 
in  different  surface  tensions  (F),  thus  influencing  the  water  trans¬ 
port.  The  values  of  surface  tension  and  contact  angle  used  in  this 
model  are  listed  in  Table  1. 


2.4.3.  Phase  change 

Liquid  water  is  considered  to  appear  in  the  fuel  cell  when  the 
water  vapor  pressure  reaches  its  saturated  value  at  the  cell  oper¬ 
ating  temperature.  Contrarily,  liquid  water  is  evaporated  when  the 
water  vapor  pressure  is  smaller  than  its  saturated  value.  Water  con¬ 
densation  and  evaporation  phase  change  is  also  an  important  factor 
to  determine  the  presence  of  liquid  water  in  multi-phase  model.  In 
the  present  study,  the  variation  of  water  mass  flow  rate  was  due  to 
condensation  (or  evaporation)  as  described  in  Table  3. 


2.5.  Boundary  conditions 

The  summary  of  the  boundary  conditions  is  listed  in  Table  4.  A 
brief  description  is  provided  below. 


2.5.1.  Inlet  of  flow  channels 

Inlet  velocities,  fuel  and  oxidant  temperatures  and  mass  concen¬ 
tration  of  species  were  set  as  given  parameters  listed  in  Table  4.  Inlet 
velocities  can  be  preliminarily  calculated  from  the  inlet  flow  rates 
based  on  the  average  current  density  (JaVg)  and  excess  coefficient  go 
that  is  defined  as 

mU2,  supply  m02,  supply  M£_. 

<^h2  =  — - >  ^02  =  — -  (15) 

'"H2 ,  consumption  ' ' l02 ,  consumption 

Then,  inlet  velocity  expression  is  given  by  the  following  equa¬ 
tions: 


Uh2  ,  in 


=  &>h2Mh2 


FiveAsurf 

2F/7H2An 


Uo2,in  =  ti>o2M 02 


fave^surf 

4Fpo2Ain 


(16) 

(17) 


Table  4 

Boundary  conditions. 


Boundary  conditions 


^  —  t/an  in 
Vh2  —  ^2 ,  in 
^H20  =  yH20,an,in 
Tin, in  =  ht2,in 
O  =  t/cat,in 
V02  —  ^02,in 
^H20  =  ^H20,cat,in 
T:at,in  —  Fair/02,in 

du0ut,an  _ q 

~  dx 

dvu^  ayH2o,a 

dx  ’  dx 

man, out  _  q 

a  dx 

ou  out,  cat 

ayH2o,c 

dx  ’  dx 

^Tcat.out  q 


=  0 


=  0 


=  0 


dx 

0s  =  0 

^=0 
dy 

0s  =  VCeii 

90m 

dy 

90s  _  n.  90s  _  n.  90m  _  r».  90m  _ 


=  0 


Locations  of  application 


Inlet  of  the  anode  flow  channel 


Inlet  of  the  cathode  flow  channel 


Outlet  of  the  anode  flow  channel 


Outlet  of  the  cathode  flow  channel 


The  anode  terminal 


The  cathode  terminal 
External  boundaries 


2.5.2.  Outlet  of  flow  channels 

The  fully  developed  condition  was  assumed  to  be  applied  for  the 
velocity  field  and  species  concentrations. 

2.5.3.  The  terminal  collectors 

There  were  two  terminal  collectors  in  which  the  anode  and  cath¬ 
ode  collectors  were  connected  to  the  external  electric  circuit.  It  was 
necessary  to  determine  the  boundary  conditions  for  phase  poten¬ 
tials  on  these  interfaces. 

2.5.4.  External  boundaries 

External  boundaries  are  defined  as  all  outside  surfaces  of  PEMFC 
except  the  terminals  (the  outside  surfaces  of  the  current  collectors 
where  the  current  enters  or  departs  from).  The  zero-current-flux 
condition  was  applied  for  the  external  boundaries  due  to  no  cur¬ 
rents  coming  or  leaving. 

For  the  sake  of  simplicity,  the  cooling  channel  was  not  included 
in  the  present  model.  However,  it  is  necessary  to  control  the  cell 
temperature.  Consequently,  the  heat  transfer  was  assumed  to  take 
place  between  the  all  outside  surfaces  and  the  ambient  environ¬ 
ment  by  forced  convection.  Note  that  heat  is  produced  when  the  fuel 
cell  operates.  For  the  case  in  which  water  finally  ends  in  vapor  form, 
this  heat  generation  rate  is  defined  as  a  product  of  the  total  current 
and  potential  loss  [31]  (shown  in  Table  3).  In  order  to  ensure  that 
the  forced  convection  dissipates  the  heat  converted  from  electricity, 
the  following  condition  must  be  satisfied: 


Qconvection  —  Qgenerated 


(18) 


then  the  convective  coefficient  h  was  chosen  as 
J(1.25 -UCeii) 

AsCTs-Too) 


(19) 


It  should  be  noticed  that  it  is  necessary  to  have  an  estimated 
value  of  surface  temperature,  Tsurf,  of  the  fuel  cell  to  calculate  the 
convective  heat  transfer  coefficient  in  Eq.  (19).  In  this  study,  it  was 
assumed  that  the  temperature  difference  between  the  fuel  cell  sur¬ 
face  and  the  ambient  is  about  50  °C  when  estimated  the  convective 
heat  transfer  coefficient. 


Table  5 

The  grid  properties  of  three  different  meshes. 
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2.6.  Solution  procedure 

The  above-coupled  set  of  governing  equations  and  relative  equa¬ 
tions  were  implemented  into/Fluent®  6.2  by  developing  our  own 
User-Defined-Functions  (UDFs)  based  on  the  general  Fluent  pack¬ 
age  that  does  not  include  the  PEMFC  module  developed  by  Fluent 
Inc.  The  developed  UDFs  are  written  in  C  language  with  about  2000 
statements.  For  the  computational  grid  used  in  this  study,  there 
were  279,500  grid  cells  in  the  computation  domain  of  the  inter- 
digitated  channel  PEMFC  employed  to  simulate  the  physical  and 
electrochemical  phenomena  in  the  fuel  cell.  The  VOF  algorithm 
was  employed  to  track  the  gas-liquid  interface.  The  solution  proce¬ 
dure  for  pressure-velocity  coupling  was  based  on  PISO  algorithm 
[33]. 


3.  Results  and  discussion 

3.1.  Validation  of  grid  dependency 

The  size  of  the  mesh  (grid  resolution)  evidently  affects  the  accu¬ 
racy  of  numerical  solution.  In  other  words,  it  was  said  that  the 
results  are  related  to  the  resolution  of  the  numerical  grid.  Therefore, 
grid  dependency  would  be  done  to  check  the  solution  at  different 
grid  sizes  and  get  a  range  of  the  grid  size  in  which  the  solutions  are 
reasonably  independent  on  the  size  of  grid.  In  order  to  validate  the 
grid  dependency  in  this  study,  a  very  fine  mesh  with  high  resolution 
was  applied  to  a  straight  channel  fuel  cell  that  has  dimensions  of 
0.001  m  x  0.001  m  for  channel  cross-section  and  0.006  m  length  of 
the  channel  and  each  cell  in  the  straight  channel  has  the  same  sizes 


3  0  001 


Z.  m 


t  -  0.0003  s 


Z.  m 


•  0  0003  s 


Fig.  2.  Computational  grids  and  comparison  of  liquid  water  tracking  among  three  cases. 
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with  0.00004  m  for  three  dimensions  (Case  1  in  Table  5).  The  com¬ 
putational  domain  of  this  mesh  contains  378,750  cells.  The  other 
two  coarser  meshes  were  also  generated  on  the  same  geometry 
of  the  straight  channel  with  the  number  of  computational  cells  of 
42,000  (Case  2  in  Table  5)  and  9750  (Case  3  in  Table  5).  The  geomet¬ 


(a)  Cathode 


rical  and  grid  properties  of  the  three  cases  are  shown  in  Table  5  and 
Fig.  2.  The  coupling  among  two-phase  flow  dynamics,  VOF  interface 
tracking  process  and  electrochemical  reaction  has  been  taken  into 
account  in  the  three  cases.  For  solving  fluid,  heat,  and  mass  trans¬ 
ports  equations  coupled  with  VOF  model  and  testing  the  possibility 


Anode 


Fig.  3.  Motion,  deformation,  detachment  and  coalescence  of  water  droplets  versus  time  in  the  cathode  and  anode:  (a)  from  t= 0.50  to  0.5004  s,  (b)  from  t  =  0.501  to  0.504  s, 
(c)  from  t  =  0.56  to  0.62  s  and  (d)  t=  0.66  to  0.72  s. 
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Cathode 


t  =  0.501  s 


Anode 


t  =  0.501  s 


Fig.  3.  ( Continued ) 
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(C) 


Cathode 


t  =  0.56  s 


Anode 


t  =  0  58  s 


t  =  0.62  s 


Fig.  3.  ( Continued ) 
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t  =  0.72  s 


Anode 


t  =  0.68  s 


t  =  0.72  s 


Fig.  3.  ( Continued ). 


of  tracking  liquid  water  inside  the  channel,  a  water  droplet  was  ini¬ 
tially  added  into  the  channel  at  the  time  t  =  Os.  The  time  step  used 
in  the  simulation  for  the  three  cases  is  10-6  s.  The  simulation  is  con¬ 
sidered  to  be  completed  when  the  water  droplet  moves  out  from 
the  channel.  In  the  three  cases,  the  numerical  results  show  that  the 


motion  and  dynamic  behavior  of  the  water  droplet  can  be  inves¬ 
tigated  while  the  mass  transport  and  electrochemical  reaction  are 
simultaneously  taking  place  in  the  fuel  cell  as  shown  in  Fig.  3.  The 
results  from  these  cases  are  quite  similar  in  terms  of  the  big  picture 
of  liquid  water  behavior.  It  also  demonstrates  that  these  meshes 
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(a)  Reference  vector:  2  m/s 


(b) 


Reference  vector:  2  m/s 


llll 


Pressure  (Pa):  1000  2000  3000  4000  5000  6000  7000  8000 

0  012- 


0  02 


Pressure  (Pa):  500  1000  1500  2000  2500  3000  3500  4000 


0.02 


Fig.  4.  Velocity  vectors  and  pressure  distributions  at  t  =  0.5003  s:  (a)  on  the  middle-plane  of  cathode  channel  (Y=0.0015  m)  and  (b)  on  the  middle-plane  of  anode  channel 
(Y=  0.00317  m). 


with  variation  of  number  of  grid  cells  (grid  resolution)  could  be 
reasonable  for  solving  the  model. 

However,  it  is  noticeable  that  the  maximum  grid  resolution 
that  is  practically  feasible  is  limited  by  the  available  computa¬ 
tional  resource  and  time,  especially  if  a  large  computational  domain 
with  complex  geometry  that  is  very  popular  in  single  fuel  cell  and 
fuel  cell  stack  simulation  is  employed.  For  this  reason,  the  coarser 
meshes  are  usually  used  if  the  solution  is  good  enough  to  resolve 
the  physical  scales  of  the  model.  Hence,  the  grid  size  used  in  Case  3 
would  be  applied  on  simulation  of  the  interdigitated  channel  mesh 
in  this  study. 

3.2.  Analysis  of  liquid  water  movement  in  the  channels  and 
porous  media 

The  effects  of  liquid  water  do  not  appear  immediately  but 
instead  require  time.  The  liquid  water  effects  will  become  observ¬ 
able  as  water  accumulates.  Based  on  the  experimental  results,  it 
usually  requires  at  least  several  minutes  of  operation  to  make  liquid 
water  effects  observable  by  naked  eyes.  In  order  to  save  the  calcu¬ 
lation  time,  an  unsteady  single-phase  model  was  simulated  from 
the  time  t= 0-0.5  s.  To  initialize  the  system,  a  series  of  liquid  water 
droplets  was  suspended  on  both  the  anode  and  cathode  channels  at 
the  time  t  =  0.50  and  0.501  s  to  further  investigate  two-phase  flow 
behavior,  especially  the  liquid  water  behavior  across  the  porous 
media,  together  with  the  electrochemical  reaction,  heat  and  mass 
transfer.  Due  to  the  structure  of  the  interdigitated  channel,  the  gas 
and  liquid  water  supply  are  forced  under  the  bipolar  plate  and  into 


the  electrode,  driving  the  water  with  it  [31  ].  Therefore,  liquid  water 
transport  and  its  effects  would  be  entirely  observed  in  the  channels 
and  porous  media  as  well  by  implementing  the  VOF  technique. 

Fig.  3  shows  the  processes  of  motion,  deformation,  detachment 
and  coalescence  of  water  droplets  versus  time  on  both  the  cathode 
and  anode  channels  and  porous  media  at  four  sequential  periods 
from  t= 0.50  to  0.72  s.  Initially,  the  droplets  in  spherical  shape  were 
added  into  the  channel  at  t=0.50  and  0.501  s  (the  positions  of  the 
initial  droplets  are  shown  in  Fig.  3a). 

3.2.1.  From  t  =  0.50  to  0.5004  s  (Fig.  3a ) 

Due  to  the  droplet  surface  tension  and  shear  stress  from 
surrounding  gas  flow,  the  droplets  were  elongated  downstream  fol¬ 
lowing  the  flow  direction  to  the  branches  as  time  progressed.  If 
the  velocity  gradients  are  large  enough,  the  surface  tension  forces 
become  unable  to  maintain  the  fluid  particles  intact,  and  thus  it 
ruptures  into  smaller  particles.  Hence,  after  hitting  the  solid  wall  in 
the  channels,  these  deformed  droplets  had  the  tendency  to  detach 
and  then  were  broken  into  tiny  deformed  ones  as  slugs  scattered 
on  the  wall  surface.  This  period  takes  place  in  a  short  time  due  to 
high  flow  velocity  in  the  channels. 

3.2.2.  From  t  =  0.501  to  0.504  s  (Fig.  3b) 

The  small  droplets  sticking  on  the  turning-wall  would  slowly 
move  forward  in  the  flow  direction  due  to  the  effect  of  wall  adhe¬ 
sion  and  surface  tension  while  the  other  droplets  from  upstream 
are  continuously  approaching  the  turn,  and  further  broken  into  the 
other  smaller  ones  in  the  same  way.  At  this  moment,  coalesce  would 
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Fig.  5.  Velocity  vectors  and  pressure  distributions  at  t= 0.62  s:  (a)  on  the  cathode  GDL/catalyst  layer  interface  (Y=  0.0023  m)  and  (b)  on  the  anode  GDL/catalyst  layer  interface 
(Y=  0.00237  m). 


be  presented  in  the  droplets,  making  a  high  concentration  of  liquid 
water  in  the  turning  area.  Consequently,  the  water  droplet  grows 
bigger.  The  shear  stress  due  to  the  increasing  air  velocity  through 
the  blocked  channel  would  increases  and  then  pushes  the  water 
droplets  downstream  and  elongates  the  droplets  into  “water  bands” 
at  the  corners.  However,  these  did  not  last  long.  Due  to  strong  shear 
stress  of  the  main  flow  in  the  straight  channel,  the  water  bands  were 
broken  again  into  the  slugs  that  would  gradually  move  to  the  end 
of  the  channels. 

3.2.3.  From  t  =  0.56  to  0.62s  (Fig.  3c) 

At  this  period,  the  slugs  gradually  moved  to  the  end  of  the  chan¬ 
nel,  then  built  up  at  the  regions  under  the  channel-ends  in  both  the 
cathode  and  anode  porous  media.  Initially,  the  water  slugs  were 
pushed  to  the  porous  media  from  the  channel  by  the  secondary 
flow  (the  flow  in  Y-direction  from  the  channel  to  the  porous  media 
applied  for  both  the  cathode  and  anode).  In  the  porous  media, 
those  slugs  were  also  strongly  influenced  by  the  primary  flow  (the 
flow  in  X-Z  planes  that  is  distributed  in  the  porous  media).  For 
slugs  in  the  porous  media  of  an  interdigitated  channel,  shear  force 
caused  by  flow  velocity  dominates  the  removal.  Depending  on  the 
flow  velocity,  the  slugs  would  be  deformed  and  spread  to  annu¬ 
lar  films  in  the  cathode  porous  media  or  would  keep  their  shape 
unchanged  as  slugs  in  the  anode  porous  media.  In  other  words, 


it  can  be  said  that  the  different  velocity  distributions  of  the  gases 
would  influence  the  formation  of  liquid  water  in  the  porous  media. 
Note  that  the  inlet  velocity  in  the  cathode  is  higher  than  that  in  the 
anode. 

3.2.4.  From  t  =  0. 66  to  0. 72  s  (Fig.  3d ) 

In  the  cathode  porous  media,  the  water  slugs  were  elongated 
by  a  strong  shear  force  as  the  time  progressed.  This  elongation  was 
continued  until  the  surface  tension  and  wall  adhesion  was  small 
than  the  shear  force  of  the  flow,  then  the  slugs  detached  to  smaller 
droplets  ( t  =  0.70  s).  Those  droplets  tend  to  be  removed  to  the  outlet 
channels  by  the  gas  flow  and  no  longer  appear  in  the  cathode  cat¬ 
alyst  layer  and  thereafter,  the  GDL.  In  fact,  liquid  water  removal  is 
feasible  in  the  porous  media  of  the  interdigitated  channel  in  which 
the  flow  field  was  forced  to  pass  through  the  porous  media,  result¬ 
ing  in  strong  convection  in  this  region.  In  other  types  of  flow  channel 
such  as  serpentine,  parallel  or  serpentine-parallel  channels,  the 
flow  distribution  in  the  porous  media  is  usually  much  weaker  when 
applying  the  same  inlet  mass  flow  rate  or  velocity.  As  a  result,  the 
liquid  water  is  difficult  to  be  removed  in  such  types  of  channel  PEM- 
FCs.  Even  in  the  interdigitated  channel  with  low  inlet  velocity  (as 
the  condition  used  for  the  anode  channel  of  this  PEMFC),  the  liq¬ 
uid  water  still  remains  unchanged  in  the  porous  media  as  shown  in 
Fig.  2c  and  d  from  t  =  0.56  to  0.72  s. 
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3.3.  Velocity  and  pressure  distributions  under  the  presence  of 
liquid  water 

3.3.1.  In  the  channels 

Fig.  4a  and  b  presents  velocity  vectors  and  pressure  distribu¬ 
tions  on  the  middle-plane  of  the  cathode  and  anode  channels  (at 
Y=  0.0015  and  0.00317  m),  respectively,  at  the  time  t= 0.5003  s.  In 
the  interdigitated  channels,  the  flow  field  shown  in  those  figures 
is  in  a  maze  with  no  outlet.  The  reactant  gases  are  driven  by  the 
channel  pressure,  forced  under  the  collector  ribs  and  into  the  elec¬ 
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trode.  It  can  be  clearly  seen  that  the  velocity  vectors  in  the  X-Z 
plane  are  highly  distributed  at  the  inlet,  decreasing  along  the  main 
channel  to  branches  and  getting  very  small  at  the  channel  cut  (the 
end  of  the  branches).  The  red  contours  indicate  the  liquid  water 
that  was  deformed  due  to  a  combination  of  droplet  surface  ten¬ 
sion  and  shear  stress  from  surrounding  gas  flow.  Therefore,  the 
position  and  shape  of  deformed  droplets  distributed  at  the  specific 
time  are  quite  different  based  on  the  inlet  velocities  at  the  cathode 
and  anode  channels.  In  other  words,  the  presence  and  blockage  of 
water  droplets  over  the  channel  led  to  the  flow  field  distributing 
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Fig.  6.  Time  evolution  of  the  motion  of  liquid  water  in  the  cathode  GDL  and  catalyst  layer:  (a)  from  t=  0.54  to  0.60  s  and  (b)  from  t  =  0.62  to  0.68  s. 
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Fig.  6.  ( Continued ). 


in  various  ways,  depending  on  where  the  droplets  were  located.  In 
addition,  pressure  distributions  in  the  interdigitated  channels  are 
insignificantly  changed  throughout  the  channel  due  to  the  fact  that 
the  flow  resistance  in  the  channel  is  very  small  compared  with  the 
resistance  in  the  porous  media.  More  interestingly,  with  the  pres¬ 
ence  of  droplets,  the  pressure  drop  reaches  its  highest  values  at 
around  the  location  of  the  droplets  due  to  the  momentum  resistance 
caused  by  the  liquid  water  in  the  channel.  It  can  be  observed  that 
the  peaks  of  pressure  coincided  with  the  locations  of  droplets  in  the 
channel. 


3.3.2.  In  the  porous  media 

Liquid  water  is  initially  added  into  the  channel,  gradually  moved 
in  the  channels  and  went  through  the  porous  media  at  the  end 
of  channel  branches  of  both  electrodes  as  shown  in  Fig.  3d.  Sim¬ 
ilar  to  those  in  the  channels,  pressure  and  velocity  distributions 
in  the  porous  media  are  dramatically  and  insignificantly  changed 
at  the  location  where  liquid  water  accumulates.  Without  the  pres¬ 
ence  of  liquid  water,  the  velocity  distributions  in  the  porous  media 
are  influenced  by  the  corresponding  position  of  the  interdigitated 
channel.  Velocity  vectors  tend  to  go  through  the  porous  media  from 
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Fig.  7.  Oxygen  and  water  vapor  concentration  distributions  on  the  cathode 
GDL/catalyst  layer  interface  (7=0.0023  m)  at  t  =  0.62  s. 
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Fig.  8.  Hydrogen  and  water  vapor  concentration  distributions  on  the  anode 
GDL/catalyst  layer  interface  (Y=  0.00237  m)  at  t  =  0.62  s. 


the  inlet  channel  to  the  outlet  channel,  such  is  a  purpose  of  the 
interdigitated-shaped  channel  that  mentioned  above.  However,  it 
would  be  different  if  there  was  liquid  water  somewhere  in  the 
porous  media.  One  can  notice  that  unlike  the  channel,  the  veloc¬ 
ity  of  the  gas  flow  in  the  porous  region  is  quite  low  (approximately 
at  the  order  of  10-1  m  s-1  or  under),  leading  to  a  small  shear  stress 
that  could  not  constrain  the  surface  tension  and  adhesion  of  liquid 
water  buildup  in  this  region.  In  other  words,  it  is  likely  that  the  liquid 
water  in  the  porous  media  barrier  of  the  gas  flow,  result  in  sudden 
changes  of  the  flow  distribution  in  the  locations  occupied  by  the 
liquid  droplets,  as  shown  in  Fig.  5a  and  b.  In  addition,  liquid  water 
also  causes  a  significant  pressure  drop  inside  the  porous  media.  In 
fact,  the  pressure  difference  is  a  driving  force  that  causes  the  reac¬ 
tant  gases  to  flow  through  the  GDLs  and  catalyst  layers  from  the 
inlet  channels  to  outlet  channels.  This  difference  would  be  much 
larger  if  the  liquid  water  occupied  somewhere  in  the  flow  field  of 
the  porous  media  blocks. 

Fig.  6a  and  b  shows  the  time  evolution  of  the  motion  of  liquid 
water  in  the  cathode  GDL  and  catalyst  layer  when  a  portion  of  liq¬ 
uid  water  emerges  from  the  porous  media.  The  gas  flow  from  the 
channel  drags  liquid  water  towards  the  end  of  the  channel  branches 
and  consequently  pushes  deformed  droplets  through  the  cathode 
GDL  under  the  channel  as  seen  in  Fig.  6a  at  t  =  0.54  s.  It  is  noticed 
that  most  of  the  liquid  water  emerges  from  the  porous  media  at  the 
channel-end  due  to  the  fact  that  there  is  no  flow  exit  of  the  inlet 
channel.  As  mentioned  above,  the  gas  flow  velocity  in  the  GDL  is 
small  due  to  porous  resistance,  hence  decreasing  shear  stress  acts 


on  the  surface  of  the  water  droplets.  Intuitively,  this  shear  stress 
seems  to  be  much  weaker  than  that  in  the  channels.  As  a  result, 
the  growth  and  spread  of  water  droplets  in  the  porous  media  has 
a  longer  time  evolution  compared  to  the  channels  as  shown  in  Fig¬ 
ures  from  t=0.54  to  0.66  s.  In  other  words,  the  growth  and  spread 
of  water  droplets  are  quite  slow  in  the  GDL  as  well  as  the  catalyst 
layer.  In  Fig.  6a  and  b,  it  also  can  be  clearly  seen  that  there  is  a 
large  amount  of  liquid  water  concentrated  in  the  region  under  the 
channel-end.  Initially,  liquid  water  is  pushed  from  the  channel  to 
the  porous  media  by  the  secondary  flow  in  the  Y-direction,  the  liq¬ 
uid  water  region  is  formed  with  a  circular  shape  due  to  uniform 
velocity  distribution.  Furthermore,  this  region  is  influenced  by  the 
primary  flow  field  distributed  in  the  porous  media  ( X-Z  planes) 
instead  of  the  secondary  flow.  The  radial  distribution  of  the  flow 
velocity  in  the  GDL  region  under  the  channel-end  would  make  the 
liquid  region  deformed  and  elongated  in  the  Z-direction  as  time 
progressed.  Meanwhile,  a  small  portion  of  liquid  water  in  the  GDL 
is  pushed  to  the  cathode  catalyst  layer  by  secondary  flow  in  the 
cathode  GDL,  and  also  is  deformed  under  the  primary  flow  in  the 
catalyst  layer.  The  presence  of  too  much  liquid  water  in  the  GDL 
and  especially  in  the  catalyst  layer  would  severely  affect  the  fuel 
cell  performance.  A  large  amount  of  liquid  water  liquid  should  be 
somehow  removed  from  the  porous  media  to  prevent  flooding  but 
a  small  amount  should  be  kept  to  maintain  a  high  conductivity  of 
the  membrane.  Interestingly,  the  simulation  results  show  that  it 
is  feasible  to  remove  the  liquid  water  by  a  high  velocity  field  in 
the  porous  media  of  an  interdigitated  channel  PEMFC  as  compared 
to  the  serpentine  or  parallel  channel  PEMFCs.  The  deformed  liq- 
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uid  water  would  continuously  elongate  until  it  is  partially  removed 
to  the  outlet  channel.  This  is  illustrated  by  a  decrease  of  the  vol¬ 
ume  fraction  of  liquid  water  in  the  liquid  region  shown  in  Fig.  6b  at 
t=  0.66  s.  It  could  be  more  clearly  seen  in  the  catalyst  layer  region 
that  a  large  amount  of  liquid  water  first  is  pushed  to  the  catalyst 
layer  from  the  GDL  at  t  =  0.54  s,  then  it  is  deformed  and  elongated 
under  the  impact  of  gas  flow  from  t= 0.56  to  0.62  s,  and  eventually 
is  taken  away  from  the  catalyst  layer  at  t= 0.64  s. 

3.4.  Species  concentration  distributions  in  the  porous  media 
under  the  presence  of  liquid  water 

3.4.1  In  the  cathode  GDL/ catalyst  layer  interface 

Fig.  7  shows  the  distributions  of  oxygen  and  water  vapor  con¬ 
centrations  on  the  cathode  GDL/catalyst  layer  interface.  The  oxygen 
concentration  is  quite  high  in  the  porous  region  under  the  inlet  and 
outlet  channels  and  the  area  between  the  channels  due  to  a  strong 
flow  convection  in  the  region.  Contrary  to  this,  the  oxygen  con¬ 
centration  in  the  area  adjacent  to  the  edges  is  lower  due  to  a  low 
convection.  It  is  realized  that  the  flow  field  in  the  area  near  the 
edges  is  insufficient  and  the  oxygen  concentration  is  dominantly 
due  to  diffusion  terms.  Under  the  presence  of  liquid  water  inside 
the  porous  media,  the  oxygen  concentration  is  significantly  influ¬ 
enced.  As  shown  in  Fig.  7,  there  are  few  spots  in  the  porous  media 
that  indicated  low  values  of  oxygen  or  water  vapor  concentrations. 
Undoubtedly,  those  spots  represent  the  regions  that  were  occupied 
by  liquid  water.  The  liquid  water  blocks  the  mass  transport  of  gas 
species,  resulting  in  a  low  concentration  on  the  distribution  of  both 
oxygen  and  water  vapor  that  were  consumed  and  generated  in  the 
cathode  catalyst  layer,  respectively.  Interestingly,  water  vapor  can 
be  seen  to  be  built-up  in  the  area  near  the  edges.  It  is  the  same  phe¬ 
nomenon  for  the  oxygen  concentration,  water  vapor  generated  by 
chemical  reaction  in  the  catalyst  layer  is  not  easy  to  remove  with 
a  strong  flow  field  in  the  area  adjacent  to  the  edges  as  discussed 
above  in  the  previous  section. 
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3.4.2.  In  the  anode  GDL/catalyst  layer  interface 

Fig.  8  shows  the  distributions  of  hydrogen  and  water  vapor  con¬ 
centration  on  the  anode  GDL/catalyst  layer  interface.  Similar  to  the 
oxygen  in  the  cathode,  hydrogen  is  also  considered  as  a  reactant 
gas  in  the  anode.  Flowever,  the  distribution  of  the  hydrogen  con¬ 
centration  in  the  interface  is  quite  different  from  the  oxygen  in  the 
cathode.  Flydrogen  seems  to  be  evenly  distributed  on  the  interface 
regardless  of  the  distribution  of  the  flow  field  (convection  effect)  in 
the  porous  media.  This  can  be  clearly  explained  by  the  high  diffu- 
sivity  of  hydrogen  gas  which  causes  the  effect  of  convection  on  the 
hydrogen  concentration  to  be  negligible.  The  way  that  liquid  water 
effects  the  hydrogen  distribution  is  similar  to  that  effects  the  oxy¬ 
gen  distribution  on  the  cathode  GDL/catalyst  layer  interface.  The 
build-up  of  water  vapor  in  the  area  adjacent  to  the  edges  also  can 
be  seen  in  the  anode  catalyst  layer  and  GDL.  It  is  noticeable  that 
the  water  vapor  is  supplied  to  the  anode  inlet  as  a  humidifier  and 
it  would  be  dragged  through  the  membrane  from  the  anode  to  the 
cathode  due  to  electro-osmotic  drag  force.  The  amount  of  dragged 
water  depends  on  the  reaction  rate  or  in  other  words  the  current 
density.  Therefore,  it  could  not  be  doubted  that  the  water  vapor  con¬ 
centration  in  the  porous  region  under  the  inlet  and  outlet  channels 
and  the  area  between  the  channels  is  lower  than  the  other  regions 
due  to  a  high  current  density  distributed  in  the  region. 

3.5.  Temperature  distribution  under  the  presence  of  liquid  water 

Fig.  9  shows  the  temperature  distributions  on  the  cathode  and 
the  anode  GDL/catalyst  layer  interfaces,  respectively.  The  heat  gen¬ 
eration  is  due  to  the  electrochemical  reaction  and  Ohmic  heating. 
Therefore,  the  temperature  is  proportional  to  the  reaction  rate, 


Fig.  9.  Temperature  distributions  at  t  =  0.62  s:  (a)  on  the  cathode  GDL/catalyst 
layer  interface  (Y=  0.0023  m)  and  (b)  on  the  anode  GDL/catalyst  layer  interface 
(Y=  0.00237  m). 

which  corresponds  to  the  local  current  density  distribution  in  the 
membrane  (Fig.  11).  Note  that  the  heat  transport  in  fuel  cells  is  by 
diffusion,  conduction  and  convection.  The  temperature  distribution 
through  the  solid  materials  (including  flow  plates/collectors  and 
solid  phase  of  porous  media)  is  also  performed  by  the  heat  conduc¬ 
tion  process.  Depending  on  different  solid  materials,  the  thermal 
conductivity  may  change,  resulting  in  different  temperature  distri¬ 
butions  in  different  solid  layers.  A  uniform  temperature  distribution 
in  the  porous  media  is  due  to  the  high  thermal  conductivity  of  the 
solid  material.  Flowever,  the  lowest  temperature  affected  by  con¬ 
vection  could  be  observed  at  the  flow  inlet  where  the  gas  flows  in 
at  the  ambient  temperature.  More  interestingly,  in  the  two-phase 
model,  the  presence  of  liquid  water  also  affects  the  distribution  of 
the  cell  temperature  by  absorbing  heat  from  the  cell,  locally  cooling 
down  the  temperature  of  the  region  where  liquid  water  occupied, 
as  seen  in  Fig.  8.  As  mentioned  above,  by  controlling  the  convective 
coefficient  h  (Eq.  (18)),  the  cell  temperature  distribution  is  in  a  good 
range  of  325-337  °C. 

3.6.  Water  content,  water  flooding  and  their  effects  on  the  current 
density  distributions  in  the  porous  media  under  the  presence  of 
liquid  water 

Fig.  10  presents  water  content  on  the  cathode  GDL/catalyst  layer 
interface.  The  water  content  is  a  function  of  water  activity  [29] 
which  is  the  ratio  of  the  water  vapor  pressure  to  the  saturated  pres- 
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Fig.  10.  Water  content  distributions  on  the  cathode  GDL/catalyst  layer  interface 
(Y=  0.0023  m)  at  t  =  0.62  s. 


sure.  In  the  porous  media,  water  content  and  ionic  conductivity 
are  strongly  related.  In  general,  the  proton  conductivity  increases 
linearly  with  increasing  water  content.  High  water  content  in  the 
porous  media  is  desired  to  increase  the  conductivity  and  the  current 
density  as  well.  Similar  to  the  distribution  of  water  vapor  concen- 


I  I  i 


Current  density  (AAn2):  1000  2000  3000  4000  5000  6000  7000 


■  i  i  1  |  i  ■  ■  i  |  ■  i  i 

0.005  0.01  0.015  0.02 

X,  m 


t  r 


(b)  Current  density  (AAn2):  0  2000  4000  6000  8000  10000 


Fig.  11.  Average  current  density  distributions  at  t= 0.62  s:  (a)  on  the  cathode 
GDL/catalyst  layer  interface  (Y=  0.0023  m)  and  (b)  on  the  mid-plane  of  the  mem¬ 
brane  (Y=  0.002335  m). 
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Fig.  12.  Liquid  water  amount  occupied  in  the  cathode  and  anode  catalyst  layers  and 
average  current  density  variations  with  respect  to  time. 


tration,  it  can  be  explained  for  the  high  water  content  in  the  porous 
media  under  the  area  adjacent  to  the  edges  and  the  low  water  con¬ 
tent  in  the  area  under  the  channel  in  which  water  is  taken  away 
by  forced  convection.  Contrary  to  this,  the  presence  of  liquid  water 
severely  directly  influences  on  the  current  density  by  blocking  the 
mass  transport  of  reactant  gas,  restricting  reactant  availability  and 
degrading  the  chemical  reactions  in  the  catalyst  layers,  resulting  in 
a  low  current  density  distribution  or  performance  loss  where  liq¬ 
uid  water  exists  as  seen  in  Fig.  11.  Intuitively,  Fig.  11  shows  that 
the  current  density  where  the  flooding  occurs  was  observed  at  a 
value  under  1000 Am-2.  Meanwhile,  an  average  current  density  is 
valued  at  6000 Am-2  and  high  values  are  over  7000 Am-2.  Fur¬ 
thermore,  it  can  be  noticed  that  the  liquid  water  distribution  in 
the  porous  media  is  not  generally  uniform  and  it  might  be  easy  to 
remove  under  a  high  convection  due  to  the  structure  of  interdigi- 
tated  channel  design.  The  removal  of  water  from  under  the  lands 
into  the  channel  is  essential  to  avoid  flooding.  However,  a  dry  con¬ 
dition  can  result  in  low  performance  and  durability.  It  could  then 
be  to  trap  liquid  water  somehow  in  order  to  humidify  the  porous 
media,  especially  in  the  anode.  Hence,  the  interdigitated  channel 
might  have  a  limitation  in  keeping  the  liquid  water  for  this  purpose. 

Fig.  12  shows  the  amount  of  liquid  water  present  in  the  cath¬ 
ode  and  anode  catalyst  layers  and  the  average  current  density  with 
respect  to  time.  The  time  evolution  is  taken  from  t= 0.52  to  0.72  s 
when  the  average  current  density  starts  to  be  stable  until  it  nearly 
becomes  constant.  In  the  anode,  liquid  water  seems  to  continuously 
penetrate  to  and  build-up  in  the  catalyst  layer  without  removal  as 
discussed  in  the  previous  section.  In  the  cathode  catalyst  layer,  this 
can  be  divided  into  three  periods:  liquid  water  penetrates  to  the  cat¬ 
alyst  layer  (t= 0.52-0.56  s),  liquid  water  accumulates  in  the  catalyst 
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layer  (t= 0.58-0.61  s)  and  liquid  water  is  gradually  removed  from 
the  catalyst  layer  (t  =  0.61-0.72  s).  Correspondingly,  the  average  cur¬ 
rent  density  of  the  cell  has  changes  due  to  the  effects  of  flooding  and 
the  catalyst  layer  humidification.  The  average  current  density  con¬ 
tinuously  increases  during  the  three  periods  although  the  flooding 
has  taken  place  in  the  cathode  catalyst  layer.  It  could  be  explained 
that  the  water  vapor  accumulated  in  the  porous  media  by  humid¬ 
ification  and  generation  in  the  cathode,  resulting  in  a  high  ionic 
conductivity  of  the  catalyst  layer  and  the  membrane,  and  increasing 
the  current  density  distribution  over  the  active  surface.  The  flood¬ 
ing,  slightly  influences  the  current  density  by  degrading  the  rising 
slope  of  current  density  (as  shown  during  the  second  period).  In 
this  case,  however,  it  seems  to  be  insignificant  in  terms  of  the  cell 
performance. 

4.  Conclusion 

In  this  study,  the  motion  of  water  droplets  in  the  channels  and 
through  the  porous  media  has  been  investigated  to  give  a  good 
understanding  of  liquid  water  behavior  including  deformation,  coa¬ 
lescence  and  detachment.  The  trace  of  liquid  water  in  the  channel 
and  porous  media  was  also  predicted  by  analyzing  the  character¬ 
istics  of  liquid  water  and  the  influence  of  the  interdigitated  flow 
field  on  the  liquid  water  distribution.  Consequently,  the  numeri¬ 
cal  results  of  this  model  show  that  the  presence  of  liquid  water 
in  the  interdigitated  channels  and  porous  media  directly  and  sig¬ 
nificantly  influences  flow  field,  species  concentrations,  heat  and 
current  transport.  The  main  conclusions  can  be  derived  as  follows: 

1.  A  comprehensive  3-D,  unsteady  PEMFC  model  combined  with 
a  VOF  interface  tracking  algorithm  is  a  powerful  and  effective 
tool  available  for  predicting  the  motion  of  liquid  water  in  the 
channel  and  porous  media  of  a  PEMFC  and  investigating  liquid 
water  effects  on  the  cell  performance. 

2.  Liquid  water  removal  is  feasible  in  the  channel  and  porous  media 
of  the  interdigitated  PEMFC.  The  reactant  flow  and  transport 
to  the  channel  and  porous  media  in  the  interdigitated  chan¬ 
nel  PEMFC  is  primarily  subject  to  forced  convection.  Hence,  this 
design  effectively  enhances  the  liquid  water  removal  and  forces 
the  reactant  flow  to  the  GDLs  and  catalyst  layers.  The  results 
show  that  the  removal  of  liquid  water  strongly  depends  on  the 
magnitude  of  flow  field. 

3.  The  presence  of  liquid  water  in  the  channel  and  porous  media 
induces  an  uneven  gas  distribution  and  creates  the  high-pressure 
regions  in  which  the  liquid  water  occupies.  As  a  result,  this  would 
provide  inadequately  distributed  mass  transfer  to  the  catalyst 
layers  and  requires  a  high  pressure  supply  for  the  fuel  cell— in 
other  words,  it  would  consume  more  power  for  fuel  cell  opera¬ 
tion. 

4.  By  investigating  the  behavior  of  liquid  water  in  the  catalyst  layers 
at  different  periods  of  time,  it  can  be  clearly  observed  that  the 
liquid  water  blocks  the  gas  transport  in  the  fuel  cell,  resulting  in 


a  degradation  of  local  current  density.  This  severely  affects  the 
cell  performance,  especially  when  the  flooding  is  significant. 

5.  The  mass  transport  of  reactant  gas  could  be  enhanced  by  forced 
convection  in  the  interdigitated  design,  however  the  water 
vapor  build-up  inside  the  porous  media  under  the  collector  ribs 
(between  the  inlet  and  outlet  channels)  can  feasibly  be  taken 
away.  Hence,  this  effect  will  reduce  ionic  conductivity  of  the 
porous  media  and  eventually  degrade  the  cell  performance.  This 
might  be  one  of  the  main  disadvantages  of  the  interdigitated 
PEMFC  design. 
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